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Abstract. Ordinary differential equations (ODEs) and ordinary difference 
systems (OASs) invariant under the actions of the Lie groups SLa;(2), SLy(2) 
and SLa;(2) x SLy(2) of projective transformations of the independent variables 
X and dependent variables y are constructed. The ODEs are continuous limits 
of the OASs, or conversely, the OASs are invariant discretizations of the ODEs. 
The invariant OASs are used to calculate numerical solutions of the invariant 
ODEs of order up to five. The solutions of the invariant numerical schemes are 
compared to numerical solutions obtained by standard Runge-Kutta methods and 
to exact solutions, when available. The invariant method performs at least as well 
as standard ones and much better in the vicinity of singularities of solutions. 


1. Introduction 

The application of Lie groups to the study of difference equations is a relatively new 
topic that has been actively pursued for the last 30 years or so. For recent reviews we 
refer the reader to [siiiiiiniiiiiEiiisQiiiaiis]. Some of the original articles pertinent 

for this study are [TIIIIIIIIIIIIIIIITIIIIIIMIEZIISH]- 

This line of research has several aspects. From the point of view of physics, 
one aim is to preserve such fundamental symmetry properties as Lorentz, Galilei and 
conformal invariance in a discrete space-time. From the point of view of mathematics, 
both pure and applied, the aim is to turn Lie group theory into an efficient tool for 
studying the solution space of difference equations as it has long been for differential 
ones H ini EH IMl- From the point of view of computing, this approach belongs 
into the field of geometrical integration [H US [35]. The aim is to improve the 
qualitative and quantitative features of numerical solutions of differential equations 
by introducing difference systems that have the same Lie point symmetry groups as 
their continuous limits (invariant discretization). 
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The original idea [7] has been applied to both ordinary differential equations 
(ODEs) [21 HU HSl mi mi mi IMI SOl and partial differential equations (PDEs) 

[iiiiiiimiiiiiiimiiiiiMimiiMiiii]- 

For first order ODEs the method provides exact discretizations, i.e. differential 
systems that have the same solutions as the ODEs m- For second and higher order 
ODEs invariant discretization often provides difference schemes that can be solved 
analytically, using either a Lagrangian [121 HI] or the adjoint equation method m- 
It was shown for numerous second and third order ODEs that invariant discretization 
provides qualitatively better fits to solutions than standard methods, specially in the 
neighborhoods of singularities [2112H| • 

The purpose of this article is to extend the method of invariant discretization 
of differential equations to larger Lie groups and higher order ODEs than have been 
treated so far. More specifically, we consider the direct product group SLa;(2) x SLy(2), 
where x and y are the independent and dependent variables, respectively and treat 
ODEs up to order 5. 

In Section 2 we briefly outline the general method of invariant discretization for 
ODEs. In Section 3, we sum up the differential invariants of SLy(2) (to all orders), 
and of SL2,(2) and SL2,(2) x SLj^(2) (up to order 5). The main results of the article 
are presented in Section 4. Thus we derive complete sets of difference invariants up 
to order 5 (using 6 points on a stencil) for all 3 groups under consideration and show 
how to obtain the differential invariants in the continuous limit. Section 5 is devoted 
to numerical examples in which we compare results using the invariant discretization 
with standard numerical methods. 


2. Differential and difference invariants of a Lie group 


Let us consider a Lie group of local point transformations acting on a Euclidean plane 
with Cartesian coordinates {x,y) generated by a Lie algebra of vector fields of the 
form 


X = f.{x,y)d^ + (j){x,y)dy (1) 

We can view the Lie group as acting either on solutions (u = f(x)) of an ODE: 

^ = 2 /^^^ - F{x,y,y',y", ■ ■ .,y^^~'^'>) = 0 ( 2 ) 


or on solutions of an ordinary difference scheme OAS 

Ea = Ea{{xk,yk}k=o) = a =1,2. (3) 

The difference scheme © that we use consists of two equations, each connecting the 
N + 1 points, satisfying 


d{E^,E2) 

d{xN,yN) ^ 


( 4 ) 


so that we can calculate {xN,yN) if the previous points {xk,yk) are known. In the 
continuous limit we put 


Xn Xji—\ — C ^ 0 (5) 

where are some constants of the order q:„ ^ 1 and we require 

{El = 0} ^ {E = 0}, {E 2 = 0} ^ {0 = 0}, e ^ 0 (6) 

Thus the lattice equation goes into an identity and the difference scheme goes into the 
target ODE. 
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We will take a given group G and find a basis for the differential invariants up to 
a certain order (i.e. fix the order N of the highest derivative) and for the difference 
invariants up to the corresponding order (i.e. fix the number of points of the lattice to 
be used as iV + 1). The ODEs and the OASs will be constructed out of the respective 
invariants, Ii,..., Ik or /f*,..., respectively and we have K < L. 

Equations written in terms of these invariant will be “strongly invariant”. Other 
equations may be “weakly invariant”, i.e. invariant on some submanifold. We will 
encounter such cases below for both ODEs and OASs. 

Basically two methods exist for calculating invariants of a given group action on 
a homogeneous manifold. One is the infinitesimal method based on the prolongations 
of the vector fields representing the Lie algebra of the group [SH |32]. The other 
method is a global one, called the method of moving frames [6l [El |33j . In the second 
method it is necessary to express the group parameters in terms of the values of a 
sufficient number of the transformed variables on some section of the generic orbits. 
For simple and semisimple groups this typically leads to algebraic equations to solve. 
In particular for the SL2,(2) x SLy(2) action studied in this article that leads to a third 
order algebraic equation. We find the infinitesimal method more convenient for the 
problem at hand and we use it throughout the article. 

The group G and the vector helds A of Q act on the variables {x, y) and on 
functions y = f{x) in the same manner, whether we are considering differential 
equations or difference systems. However in the continuous case we prolong to actions 
on derivatives in a a standard manner |32j . In the discrete case we write X at some 
point Xk of the one dimensional lattice and then sum over all points involved in the 
OAS 0130]: 

Pr Xn = 'y ^ {f,{i^n+kTyn+k)dx^j^^ + (t>{Xn+k, yn+k)dy„^k') C^) 

k 

The summation over k is over all points on one stencil. The index n labels the position 
of the stencil used in the calculation. 

In both cases we find the invariant by solving the system of determining equations 
following from the invariance condition 


pTX^{x,y,y',...,y^^'>) =0 

(8) 

X„^{xn+k,yn+k) = 0 

(9) 


There will be more functionally independent difference invariants than differential 
ones. We will divide the difference ones into two sets; those that go into differential 
ones in the continuous limit and those that vanish in this limit. 

The connection between difference and differential invariants is established by 
using Taylor expansions of the discrete quantities. We restrict ourselves to a single 
stencil, i.e. points {xo,yo), ■ • ■, ixN,yN), and choose a point about which to develop, 
say xq. All other points are expressed as: 

k 

Xk = Xq hi, Uk = y{xk), 1 <k < N 

1^1 

and we expand all discrete invariants using the truncated Taylor series: 


( 10 ) 
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The continuous limit is taken as in ([^. The result will be expressed in terms of the 
basis of differential invariants. In general the limit may depend on the constants a„ 
in (§. These will be specific numbers once the lattice is chosen. Detailed examples 
will be given in Section 4. 


3. Differential invariants under SLy(2), SL 2 ,( 2 ) and SL 2 ,( 2 ) x SLy(2) 

In this section we restrict the group G to be SLy(2), SLa;(2) and SL 2 ;( 2 ) x SL 2 :( 2 ) 
respectively, and will present bases for all differential invariants up to order 5, though 
it would be easy to proceed to higher orders. The reason for this choice is that N = 5 
is the lowest order at which an SL 2 ,( 2 ) x SLy(2) invariant exists. 


3.1. Invariants o/SLy(2) 
ra of t 

prolongations: 


The Lie algebra of the group SLy(2) is generated by vector fields dy, ydy, y'^dy with 


pr(^)dj, = dy 


N 


N 


( 12 ) 


pr^^^y^dy = y^dy + , 


Solving the corresponding PDEs 
invariant 


y‘ 


3 fy 


n,' 0 ( o,' 


k=l 

1 ) for TV = 5 we find the lowest order differential 


(13) 


The third and fourth invariants can also be calculated directly. 

Alternatively, since the variable x is invariant we can start from the lowest order 
invariant involving derivatives of y, namely J 3 and generate a different basis of SLy(2) 
invariants, by invariant differentiation: 


dx'^ 


(14) 


All SLj,(2) differential invariants of order up to N will be functions of 

{x, T3, t/4, t/5,..., Ttv} (^ 3 ) 

We mention that J 3 is the well known Schwarzian derivative with many interesting 
applications [35j . We will use in this work the first three invariants, J 3 , J 4 and J 5 . 
The explicit form of J 4 and J 5 are: 

,(4) y/'y/ll /I3 

_ -U .ti 

,,/3 


J4^4 = ^-4^+3^ 

y y y 


Jr, = J^' = 


y' 


^ y/2 ' yl3 yl2 ^ yl4: 


We could also use a simplified fifth order invariant, adding a multiple of J| 


J5 = J5 + 4 Jo — 


(5) 


y 
y' 


-5 


y"yO) 


ya 2 y,n 


(16) 

(17) 

(18) 
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3.2. Differential invariants o/SLa;(2) 

The Lie algebra of SLa;(2) has the basis dx: xdx, x^dx with prolongations 
pr(^)9, = dx 


N 


=xdx-^ ky^^^dy(k) 


k^l 

N 


(19) 


= x‘^dx - ^k{{k - l)y^^ + 2xy^^^)dy(k) 




It is of course equivalent to SLj,(2) and the two are transformed into each other by 
a hodograph transformation. We treat SL 3 :( 2 ) separately here since we are interested 
mainly in the action of SLa;(2) x SLy(2). 

The three lowest order invariants are (higher order ones are also computed in a 
straightforward way 0 ) 

1 / ll'" R ^ 

K3 = 


1 

(j/O 


y 


y 


y 


y 


y 


(4) 


//3 


^ a f,y"'y” , f.y' 

Ka = —^-6- Z -h 0— 

yf4 y/5 y/6 


( 20 ) 


K. = 


y 


(5) 


///2 


y' 

- 10 ^—^-4 — 

yl5 yl6 yf6 


42 


„/7 


63 y”^ 

Yy^ 


and, of course, y itself is an invariant. 

It is interesting to consider the behavior of the SLy(2) invariant J 3 . We have 

(pr^^^Sj;) J3 = 0, (pr^^^iSa;) J3 = -2J3, ipr''^Y'^dx)J3 = -4a;J3(2I) 

Thus, while J 3 is not invariant under SLa;(2), the equation 

J3 = 0 (22) 


determines an invariant manifold and the equation ( 22 ) is “weakly invariant” under 
the entire group SL 2 ;( 2 ) x SLy(2). The same of course holds for K 3 = 0. 


3.3. Differential invariants ofSljx{2) x SLj,(2) 

We start from the invariants of SLj,(2), /(J3, J4, J5, ) and require that this function be 


annihilated by the vector fields (19). We find that there are no differential invariants 


of order TV < 5 and only one of order 5, namely 


J 5 5J; 


K. 






4J| 


Kl 


hKl 

AKl 


i 2 y"'y' - 2 ,y"‘^Y 


2 y'^{ 2 y'y"' - Zy'^y^'^) + 20 y' 3 y"y'"y( 4 ) _ ^y'\yD) f 


16y'3y'"3 + l2y’Y”Y'"^ - Wy”V + 9y"® 


(23) 










Symmetry preserving discretization of ODE 


6 


4. Difference invariants and their continuous limits for the groups SL2;(2), 
SLy(2) and SL^(2) x SLy(2) 


4-1. General comments. The cross-ratios 


The Lie group actions that we are considering in this article are the standard projective 
action of SL(2) on a real or complex line (the action of the Mobius group). The 
fundamental invariants of this action are well-known (and can easily be reobtained 
using the prescription §)• The lowest order invariants involve 4 points and are the 
cross-ratios (anharmonic ratios): 


Rk+3 = 
Sk+^ = 


ivk+i — yk+l){yk +2 — Vk) 
{Vk+S — yk+2){Vk+l — Vk) 
^fc+3 ^k+l){,'^k+2 ^fc) 

(:^fc+3 ^fc+2)(^fc-t-l ^fc) 


(24) 

(25) 


for SLy(2) and SLa,(2) respectively. 

For SLa;(2) X SLy(2) these are the only four-points invariants and all higher order 
invariants can be formed by shifting the four points to the right, forming e.g. the 
cross-ratios Rk+iiRk+b • • ■, and taking linear combinations of Rk+z, Rk+A ■ • ■, etc. 
The problem is to form the linear combinations that will have the chosen differential 
invariants of Section |3] as a continuous limit. 

For SL2,(2) and SLj,(2) we have further difference invariants, namely the 
dependent variables ?/„ and the independent variables Xn for SLa;(2) and SLy(2), 
respectively 

The cross-ratios Sj will be used to write invariant lattices, e.g. Sk+z = ^ or 
5^+3 = ASk+A, where A is a constant. 

The cross-ratios Rj will be expanded into power series, and using equations (101 
and (11) we relate them to differential invariants. 

Let us consider SLy(2), SLa;(2) and SLa;(2) x SLy(2) separately. 


4-2. The group SLy(2) 


Taking 4 adjacent points, Xo,Xi,X 2 ,X 2 (we simplify the notation writing Xn+k 
and expanding around xq we obtain: 


.Rs — S^ 


1 — -h 2 (hi -h hk +2 + hk+s) ( J 3 7 ( 3/11 -I- ‘2hk+2 + /i3)'^4 

-f — [ 3(6/1^ -|- 3/I2 -|- /I3 -f 8/11/12 + 4/11/13 -|- 3 /l 2 /l 3)>/5 


— I0{7hf Ahl -\-h\-\- IO/ 11/12 + 5 / 11/13 -I- 4112 / 13)^^3 

+0(/i®) 


= Xk) 


(26) 


We can make a similar expansion for two more sets of four points, Xi,X 2 ,X 3 ,Xi 
and X 2 ,X 3 ,Xi,X 3 (always expanding around xg). We see that to lowest order in h the 
difference invariant that has the correct continuous limit is 


2^(3) ^ ^_ L _ R±+3 \ . 

{xk +2 - Xk+l){xkA-3 - Xk) \ Sk+sJ' 


(27) 
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indeed for all values of k (we will only need fc = 0 , 1 , 2 ) we have 
lim = J 3 . 


In view of (26) we can define a new set of difference invariants: 
4 


= 

which can be expanded in h as: 


d3) _ A 3 ) \ 

^k+4 ^k+3 I 


(28) 


(29) 


( 4 ) 


1 


15(/ii + /l2 + ^3 + hi) 


3 4/if + 3/12 + 2/13 + hi 


-\-7h1h2 H” 6 /li/i 3 “h 5 /ii/l 4 -|- 5/12^3 H“ 4/12^4 “ 1 “ 3 /i 3 /l 4 ^ t/5 

—2 (hj — 3^2 + 3^3 — hi — 2/11/12 — hihz + 112^4 + 2/13/14) j| 

+0(/i^) (30) 

and a similar expression for ^ 5 ^^. In the continuous limit, we have, for all k (in 
particular, for the only two values we need, fc = 0 , 1 ): 


lim = J4 


(31) 


Similarly, to obtain the fifth order differential invariant in the continuous limit 
we form another set of difference invariants, namely 

5 


and 

r( 5 ) _ T , 10 

4^5 “ 


r( 5 ) _ 
4^5 - 


1 


X 5 - Xo 


(A - A) 


(32) 


hi 


5 /ii + /12 4“ fcs 4“ ^4 + fcs 

{hi + h 2 ){h 2 4- / 13 ) 


{hi 4 “ /12 4 “ /13 4 “ hi 4 “ h^){hi 4 - /12 4 “ /13 4 - hi) 
(/i2 4 - h^){h'i 4 - hi) 

{hi 4 “ /12 4 “ /13 4 “ /14 4 “ h^){h2 4 - /13 4 - /14 + /15) 


ji + o{h). 


(33) 


From (26) (e.g. for fc = 0, 1 and 2) we see that the limit of is not exactly J 5 of 

(34) 


(141 but rather a combination of J 5 and Jg 


lim lP = J 5 + IFo ^3 

hj^O 


with 


and 


VFo = lim W 

hj—^O 


(35) 


^ _ W /I _ Xj - X 3 _ {X 3 - Xi){x 2 - acp) ( 3:4 - a: 2 )(a :3 - a;i) \ 

3 \5 3:5 - ( 3:5 - a;o)(3^4 - 3 : 0 ) ( 3:5 - a:o)(3:5 - 3i)/ 

The coefficient W depends on the specific form of the lattice, and IFq is a finite number 
(not necessarily zero). 

For instance, let us consider an SL3;(2) x SLj,(2) invariant lattice given by 
S,=K, Wj 


(37) 
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where K is a, constant. This equation (371 was solved in m for arbitrary values of 
K. The result is particularly simple for K = 4. A particular solution, not contained 
in the general one is 

Xm = Am + B (38) 

and the general one is 

1 




(39) 


{A, B, C are constants). Solution (381 corresponds to a uniform lattice with hm = A 
for all m. For ([M| and ((Ml) we have 


ITo = 0, ITo = 


2A^{iA^ - bAB - B^) 

{A + B)(2A + B){3A + B){4A + B) 


(40) 


respectively. 

Thus, Wo is a definite number (and can be set equal to zero in the case (39) by 
choosing B = i(A(—5 ± -s/S?))- 

We see that it is not difficult to construct SLy(2) difference invariants of arbitrary 
orders. The challenge is to hnd an appropriate basis for these invariants that in the 
continuous limit reproduces the chosen basis of difference invariant. 

The connection between the difference and differential SLy(2) invariants of order 


3, 4 and 5 is given by equations (27), (29) and (32), respectively. 


4-3. The group SL2,(2) 

The approach used for SLy(2) must be modified since the differences hk are not 


invariant under SLa;(2). Instead of (27) we expand the SL2:(2) invariants: 

Rk+3\ 


— 

JU I Q - 


6 


ivk+s - yk)iyk+2 - Vk+i) 


1 - 


Sk +3 ) 


(41) 


and obtain 

ikfg ^ + — (3hi + 2/i2 + h^)y'K4^ 

1 


120 ' ^^2 + ^3 + + 4 /iihfc +3 + 3/12/13)?/' 

-4(19/i^ + 12 / 1 ^ + Ahl + 27/ii/?2 + llhihs + I2h2h3)y''^ K'^ 

+ 15(3/?? + 2hl + hl + 4/?i/?2 + 2/?i/? 3 + 2/?2/?3)?/"iF4 ) + 0{h^) (42) 


and similar expressions for k = 1,2. Thus we have 

= Vfc 


In analogy, we define 


^4^4 = 


yk+A — yk 




'fc+4 


-(3) \ 
k+ 3 ) > 


V/c 


(43) 


(44) 


lim = Ki, Vfc 


satisfying 


(45) 
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Finally, to obtain a fifth order differential invariant in the continuous limit we 


define 



2/5 - 2/0 

satisfying 

lim = Ki 

O-s-o 

where 





= __ 

6 \ 5 + /i 2 “h ^3 H“ ^4 H“ ^5 


(46) 


(47) 


(48) 


As in the case of SL 3 ;( 2 ) we do not obtain in the limit but have an additional 
term involving a scalar multiple, Wx,o, of a lower order invariant K^. The number 
Wx,o can be evaluated on any lattice and is 

= 2 (49) 

on a uniform lattice. 


4.4- The group SL 2 ;( 2 ) x SLy(2) 

For the group SL 3 ,( 2 ) x SLj,(2) all invariants must be constructed out of the cross¬ 


ratios Rk+ 3 , Sk+ 3 - To obtain the lowest order differential invariant given by (23) 


we need three values of fc, e.g. k = 0 , 1 , 2 . As always, the difficulty is to identify 
the combination or combinations of them that go to in the continuous limit. One 
way to do this is to expand Qk +3 = 1 — for A: = 0 , 1 , 2 into power series, and 
eliminate terms of order h'^ and using Sk -3 whenever possible. Another way is to 


inspire oneself by the continuous limit (23) and build up an invariant with the correct 


limit using discretized versions of J 5 , J 4 , J 3 and K^, K 4 K^. Both methods are quite 
laborious, even using computer algebra, and lead to the result 


(4)04)' 


l_|(5) _ (a^5 - a^4)(a;i - xq) 1 / _ 5 L\ 'L 


{Xi - X3){X2 - Xi) 7 ^( 3 ) 

or explicitly in terms of the invariants Ri and Si 
H(5) = 10 54 


L 


(3) 


4 T T 

-^3 ^4 . 


(50) 


3 ( 53(1 - 54 ) + 54)(54(1 - 55 ) + 55)(54(1 - 53)(1 - ^ 5 ) - 53 ^ 5 ) 


54(1-55)^+54(1-53)^ 


-(1 - 54)(54(1 - 53)(1 - 55 ) - 5355 ) 

V4 


-54(1- 53)(1- 55) 


Qi \ 


Q3Q3 / 


where 


Q, = l- 


Ri 

R.'' 


i = 3,4,5 


(51) 


(52) 


The continuous limit is given by (the quantity W was defined in (36)): 
lim = ( lim 

hj^O \hj ^0 /i4/l2 J 


(iJs + lim W) 


(53) 
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For the case of a uniform lattice {Sj = 4, hj = hj+i) we have IF = 0 and 
hm H(5) = = 4 - S 


h->-0 




4J| 


(54) 


with 

I ifS) 16i?5 + i?4(3i?4 + i?5 — 32) + — 5i?5 + 16) 

2(i?3-4)(i?4-4)(i?5-4) ^ ^ 

For other invariant lattices we have W = constant ^ 0. 

The overall conclusion from this section is that we have constructed all difference 
and differential invariants up to order 5 for all groups considered. We have also shown 
how to proceed to higher orders, using the fundamental invariants, specifically the 
cross-ratios. 

We now proceed to test the invariant O ASs presented in this section and numerical 
schemes for specific equations. So far, only the group SLy(2) has been used in this 
manner and only for ODEs of order 2 or 3 |S1 [5S] ■ 


5. Numerical examples 

We present in this section some representative examples of differential equations which 
are invariant under the groups SLy(2), SLa;(2) and SL2,(2) x SLy(2). Some particular 
solutions are numerically computed using the invariant discretizations studied in this 
work. In some cases the discretized scheme provides the exact solution. In other 
cases, the roundoff errors will forbid to find the solutions at some points. Increasing 
the working precision would allow to find more approximate solutions, but we cannot 
keep the same working precision beyond a certain point. These quantitative aspects 
of the theory have been studied elsewhere [3] and will not be discussed in detail in 
this work. Since our purpose is to provide some qualitative remarks on the invariant 
method, the standard numerical approach has been carried out using the software 
Mathematica, which allows a high level performance in all the examples we will discuss 
below, while being very simple to use. In each case, the program chooses the most 
appropriate method (for instance, an Adams predictor-corrector or an explicit Runge- 
Kutta method) and, whenever necessary, a variable step, looking for higher accuracy. 
In both approaches (standard and invariant) some control parameters, like working 
precision and accuracy, have also been chosen in order to get the best results. 


5.1. Example 1: fourth order equation, invariant under SLy(2) 


The differential equation: 


y 


(4) 


,//3 


» w .. 

4 tj tj I n 

+ 3^73 = cos 


(56) 

y'^ y-o 

is invariant under SLy(2), but not under any transformation of SL2;(2). Note that the 
equation can be written as: 

Jg = coscc (57) 

and a first integral is: 

Jg = sin a: -I- A (58) 

However, to check the invariant method, we will use the invariant scheme for fourth 
order differential equations (311. 
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We choose as initial data: 

y(l) = 1.0, y'(l) = -1.0, y"(l) = -2.5, ?/"'(l) = 5.0 (59) 

Using a standard numerical solution in the interval [1, 2.5], we obtain the graphics 
of Figure 



Figure 1. Example 1. Solution of equation ||56|l. Standard numerical method. 


We use a uniform lattice with step h (which like any lattice depending only on 
the variables Xi is invariant under SLy(2)) on a hve-point stencil. The equations are 

-1 


= 1 - 


Vn 


X yn-2 


Vn—Z yn—2 f {yn — 1 ^n—3)(2/n—2 2/n— 4 ) 


Un—S Vn—l xiUn — l 1/n—2)(2/n—3 2/n— 4 ) 

yn—l{yn—3 yn — 2 ) f {yn—1 yn—3){yn—2 yn—4) 

yn—3 yn—1 \{yn—l yn—2){yn—3 yn—4) 


— 2h? cos{h{n — 2 ) + ccq ) 


—2h^ cos(a:o + h{n — 2)))) 
Xn = Xq + rih 


with initial conditions taken from the numerical solution computed with a standard 
procedure. We compute the cosine function in the middle point of the stencil and 
consider three cases, with three different steps h. 

Table compares the values at several points using standard numerical methods 
(explicit Runge-Kutta) and the invariant approach (Inv) with different steps. The 
accuracy of the method improves when h goes to 0, (although the increasing number 
of steps will correspondingly increase the roundoff errors. We have used the same 
working precision in all cases). 


X 

Standard 

Inv h = 0.1 

Inv h = 0.01 

Inv h = 0.001 

1.5 

0.451089 

0.451095 

0.451084 

0.451089 

2.0 

0.310434 

0.310466 

0.310435 

0.310434 

2.5 

0.298849 

0.298885 

0.298850 

0.298849 


Table 1. Equation l |56[ |. Values of the solution at points 1.5, 2.0 and 2.5. The 
columns are the values obtained with a standard numerical procedure and with 
the invariant one with different steps, respectively. 


(60) 


To test the precision of the results we can use the distance, as mean square 
averages, between the numerical solutions computed by the invariant scheme and the 
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X = 


lEniylf''-yny 


Unyl 


(61) 


and get the results in Table at different values of the step, in the interval [0,2.5]. 
The table shows the deviation of the invariant solution from the standard one. The 
precision of the method is improved when h diminishes, as it was expected. In other 
cases whe n w e know the exact solution, we could compare with the exact solution 

in ra. 


h 0.1 0.01 0.001 

4.6729. X 10-^ 1.7866 x 10"® 7.2291 x 10"*^ 

Table 2. Equation | |56| l. Values of x- Invariant approach with different steps 
versus standard numerical approach 


In Figure the dots are the values of the numerical solutions computed with an 
invariant approach, for h = 0.001. The solid curve is the standard numerical approach. 



Figure 2. Equation ( |56[ |. Invariant approach versus standard approach for step 
h = 0.001 


5.2. Example 2: third order equation invariant under SLa;(2) 


Let us consider the differential equation: 

2/'2 y/ 2y'2 J 


c^O 


(62) 


where c is a constant. This equation is invariant under SL2,(2) (and under translations 
in y). The order can be easily reduced by one, and the general solution in terms of 
elementary functions can be written as 


y{x) = C3 + Y -arctanh {cix + C 2 ) 

A 2-parameter family of particular solutions not contained in (62) is 

y = a+ log(x - b) 
v 2c 


(63) 


(64) 
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When c = 1/2, a particular solution is: 

y(x) = log |x| (65) 

which has a singularity at x = 0. The initial values for the solution in the region x < 0 
can be chosen as 


y(-l) = 0, y'(-l) = -l, y"(-l) = -l (66) 

and the solution can be numerically computed by any standard method which will 
obviously stop at t = 0. 

We shall now apply the method of invariant discretization to calculate solutions 
of the type (631 and (65) numerically. Let us start with solution (651. The invariant 
approach is obtained with the difference and x-lattice equations (we choose a uniform 
lattice, yielding 5 = 4) 


6 


1 - 


(yn—3 yn—l)(yn—2 Vn) 

4(y„-3 - yn- 2 )(yn-l - Vn) 


(yn -1 - yn- 2 ){yn - Vn-s) 

Xn = Xo + nh 

The initial conditions are taken from the exact solution: 


1 

2 ’ 


(67) 


Xo = —1, h = 0.0001, ?/o = 0., yi = log(0.9999), ?/2 = log(0.9998) 

However, in contrast with Example 1, the difference equation is not linear in j/„, 
but quadratic. Apart from technical difficulties, one obtains two possible solutions. 
One of them provides the approximate solution. Note that the difference equation 
could be, in principle, computed beyond the stop point x = 0 but becomes complex 
if x„ is greater than 0. See Figure for a graphics of the approximate solution (dots) 
versus the exact one (solid line). 



Figure 3. Differential equation ( |62[ l x < 0. Invariant approach, h = 0.0001, 
(dots) versus exact solution (solid line). 


We can also compute the solution in the region x > 0. Starting at xq = 1 and 
using a negative step h = —0.0001, we get for the invariant numerical solution the 
graphics in Figure As in the previous case, the values become complex when 
x„ becomes negative. Then, we have reproduce in these graphics, and the two 
regions of the logarithmic solution. 

Let us now choose a different solution of (62) (with c = 2), namely 

y(x) = arctanhx (68) 

which exists only in the interval (—1,1). The initial values for this solution are: 

y(0)=0, y'(0) = l, y"(0)=0 (69) 
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Figure 4. Differential equation ( |62[ l, x > 0. Invariant approach, h = —0.0001, 
(dots) versus exact solution (solid line). 


and can be numerically computed by any standard method which will stop at a: = 1. 
The invariant approach is obtained with the same difference equation as above (67), 
but the initial conditions (which are again taken from the exact solution) are: 

Xo = —0.9, h = 0.01, yo = arctan(—0.9), yi = arctan(—0.89), y 2 = arctan(—0.88) 


As above, one of the two solutions for provides the approximate solution. See 
Figure for a plot of the invariant numerical solution (dots) versus the exact one 
(solid line). The x estimator in the interval [—0.9, 0.9] can be computed against the 



Figure 5. Differential equation l |62[ l. Invariant approach, h = 0.001, (dots) 
versus exact solution (solid line). 


exact solution, for different values of the step, see Table 

h 0.1 0.01 0.001 

X 0.145901 0.007028 0.001131 

Table 3. Equation l|62[|. x for different values of the step h. 


5.3. Example 3: third order equation invariant under SLa;(2) 


The differential equation: 


1 

1/2 


y 


32 /" 


= y 


y' 22/'2 


(70) 
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is invariant under SL 3 ;( 2 ). Apparently, the general solution cannot be constructed in 
terms of elementary functions. 

With the initial values 

2/(0) = 10, 2/'(0) = -l, 2/"(0) = -10 (71) 

the standard numerical methods stop at a; « 0.14. A singularity (or a point with 
infinite derivative) is expected. 

The invariant approach is obtained with the difference and a;-lattice equations (we 
choose a uniform lattice, yielding S' = 4). The difference equation is a polynomial of 
third degree in //« and the search for real solutions becomes rather involved. However, 
it can be done and the results appear in the tables and graphics we present. The role 
of the other branches, which can contain complex values, is not well understood. 

If the initial conditions are taken from the approximate solution: 

xo = 0., /i = 0.001, yo = 10., 2/1 = 9.9989, 2/2 = 9.9979 (72) 

the graphics of the invariant approach is given in Figure Although the difference 
equation could be solved beyond the point x = 1.4, 2 /n becomes complex (for the 
chosen branch). 



Figure 6. Invariant approach (dots) versus standard approach (solid line) for 
equation ( |70[ l. 


5.4- Example 4- invariant equation under SLa;(2) x SLy(2); a discrete exact solution 

We will consider differential equations invariant under the direct product group 
SL,(2) X SLy(2). 

The equations are of the form 

H 5 = C, y^^^ = pr( 2 yyVZ 17 / 2 ) f + 2(c + 4)y'^y'"^ 

-|(c+0 (42/'22/"V"^-62/'2/"V" + 32/"®)) (73) 

where c is a constant. The equation can be easily solved, although the general solution 
can adopt several equivalent forms. We will consider the particular solution: 

y{x) = 


1 - e^’ 


c = 0 


(74) 
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The values corresponding to this solution can be easily obtained by a standard 
numerical method with initial conditions: 

^ , fc = 0,l,2,3,4 (75) 

X — — 1 




1 — e^ 

until a:: = 0 where the function has a singularity (see Figure]^ black curve). 



Figure 7. The solution | |74[ l of equation | |73^ obtained by a standard numerical 
method. Grey line, the exact solution. 


The invariant approach is obtained with the difference and ai-lattice equations 
(we choose a uniform lattice, yielding 5 = 4). 

3i?4 + (i ?5 — 32)i?4 + 16i?5 + i?3(i?4 — 5i?5 + 16) = 0, Xn = xq + nh (76) 

where 

_ {yn—3 yn—h){yn—2 yn—A) 

iVn—A yn—b){,yn—2 yn—s) 

_ {yn — 1 yn—s){yn ?/n— 2 ) 

{yn-2 - yn-3)iyn “ 2/n-l) 

The difference equation is linear in t/„. The initial conditions are taken from the exact 
solution: 

= h = 0.1, = , A: = 0,1,2,3,4 (78) 

1 — e *= 

and the invariant solution is represented in Figure]^ 


i?4 = 


(jJn—2 yn—A)(jJn—l yn—s) 
(j/n-3 - yn-4){yn-l “ yn-2) ’ 


(77) 



Figure 8. The invariant approach (dots) versus the exact solution \7A\ of 
equation (|73[l 














Symmetry preserving discretization of ODE 


17 


It turns out that the solution provided by the invariant approach is an exact 
solution: 

Vn — Z : — 1 rih (^1^) 

1 — e^" 

This is easy to check. Compute for any four consecutive points, for instance R 3 , 
when pn = ■ We get (for a uniform lattice, a;„ = — 1 + hn) 

i?3 = 2 + e'* + e"^ (80) 

and the same expression for R 4 and R^. Substituting in we get zero. This can be 
also observed from another point of view. The expression for is a solution of the 
equation: 

Rk = a (81) 

and this provides a solution of the difference equation we are considering. 

This observation allows to study the problem in the opposite direction. Since 
Rk = a is a solution of — 0, the equation Ri = a should provide solutions of the 
difference equation and approximation of solutions of the differential equation. These 
exact solutions have been computed in [11] (for the case under study, see equation 
(5.24) of this reference). 


5.5. Example 5: invariant equation under SL2;(2) x SLy(2) 


Our final example is a discussion of the equation 


o 

11 

(82) 

and the particular solution: 

/ N 1 

y{x) = tan — 

X 

(83) 


The numerical values of this solution can be easily obtained by a standard 
numerical method, although it cannot be prolonged beyond the singularities (the first 
one greater than x = 0.1 is located at a: = 2/(57r)). We take an initial condition at 
Xq = 0.1 with the values of the function and its derivatives computed using the exact 
solution (see Figure]^: 

/ 2. \ 

y(a:o) = tanlO, (xq) = (tan - I , fc=l,...,4 (84) 

\ ^/a;=0.1 

The invariant approach is obtained with the difference and a:-lattice equations 
(we choose a uniform lattice, yielding 5 = 4): 

3i?| + (i ?5 — 32)i?4 + 16i?5 + i?3(i?4 — 5i?5 + 16) = 0, a;„ = a;o + nh (85) 


and is linear in y„. The initial conditions are also taken from the exact solution: 

1 


Xq = 0.1, pk = tan 


fc = 0,l,2,3,4 


( 86 ) 


^ 0.1 + kh^ 

The scheme is very sensitive to the step size and the fixed working precision. 
However, it is possible, using the invariant approach, to go beyond the singularities of 
the solution. We will just present a qualitative summary of results. 


The graphics in Figures 10 11 and 12 represent the solution and the invariant 


discretization for h — 0.01, h = 0.005, and h = 0.001, respectively. 

Smaller steps provide a better approximation, although the roundoff errors 
prevent us (even at the cost of greater working precisions) to go beyond a certain 
point. 
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Figure 9. Standard numerical approach versus the exact solution ( |83| l of equation 
(|82[l. The grey line corresponds to the exact solution. 



Figure 10. Invariant approximation for solution l |82[ l of equation | |83^ . Step 
h = 0.01 



Figure 11. Invariant approximation for solution l |82[ | of equation | |83| l. Step 
h = 0.005 


6. Conclusions 


The main theoretical results of this paper are contained in Section 4. We have shown 
that starting from the four-point difference invariants Ri (24) and Si (25) of SLy(2) 
and SLa;(2) we can construct difference invariants of arbitrary order, for the 3 groups 
considered in this article. In the continuous limit they approach the corresponding 
differential invariants. Explicitly we go up to order N = 5 and this provides invariant 
schemes for solving invariant ODEs of order up to five numerically. 

















Symmetry preserving discretization of ODE 


19 



Figure 12. Invariant approximation for solution l |82[ | of equation | |83| l. Step 
h = 0.001 


The numerical results are presented in Section 5. We consider several ODEs of 
order = 3,4 and 5. The main features that emerge are the following: 

(i) Invariant numerical methods and standard methods provide very similar results 
for smooth solutions. 

(ii) For solutions with singularities invariant methods provide significantly better 
results, specially close to singularities and beyond them. 

(iii) In some cases the “invariant numerical” solutions are exact (see example 4 for a 
fifth order ODE. This is always true for first order ODEs [40] and was already 
observed for some second order ones O [13] . 

Since symmetries are an essential part of any physical problem preserving them in a 
discretization is important in itself. This is true independently of whether invariant 
discretization improves numerical results. 

An open question which merits further study is that of identifying equations and 
initial or boundary conditions for which invariant methods provide exact solutions. 
Work in this direction is currently in progress. Another line of research is related 
to the study of the several solutions arising from nonlinear discrete schemes (implicit 
schemes, where the highest point pn is not defined as a unique function of the previous 
points in the stencil). 
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